Operational stability study of lactate biosensors: modeling, parameter identification, and stability analysis

Introduction This paper investigates the operational stability of lactate biosensors, crucial devices in various biomedical and biotechnological applications. We detail the construction of an amperometric transducer tailored for lactate measurement and outline the experimental setup used for empirical validation. Methods The modeling framework incorporates Brown and Michaelis–Menten kinetics, integrating both distributed and discrete delays to capture the intricate dynamics of lactate sensing. To ascertain model parameters, we propose a nonlinear optimization method, leveraging initial approximations from the Brown model’s delay values for the subsequent model with discrete delays. Results Stability analysis forms a cornerstone of our investigation, centering on linearization around equilibrium states and scrutinizing the real parts of quasi-polynomials. Notably, our findings reveal that the discrete delay model manifests marginal stability, occupying a delicate balance between asymptotic stability and instability. We introduce criteria for verifying marginal stability based on characteristic quasi-polynomial roots, offering practical insights into system behavior. Discussion Qalitative examination of the model elucidates the influence of delay on dynamic behavior. We observe a transition from stable focus to limit cycle and period-doubling phenomena with increasing delay values, as evidenced by phase plots and bifurcation diagrams employing Poincaré sections. Additionally, we identify limitations in model applicability, notably the loss of solution positivity with growing delays, underscoring the necessity for cautious interpretation when employing delayed exponential function formulations. This comprehensive study provides valuable insights into the design and operational characteristics of lactate biosensors, offering a robust framework for understanding and optimizing their performance in diverse settings.


Introduction
1.1 Operational stability of biosensor vs. Lyapunov stability of the dynamic model Operational stability of biosensors means "retention of activity of a protein or enzyme when in use" (Gibson, 1999).
It corresponds mainly with the same notion for dynamic systems from stability theory.Traditionally, enzyme-substrate interaction is simulated with the help of the Michaelis-Menten model, which is a nonlinear dynamic system.On the other hand, previous studies have paid little attention to the qualitative behavior of the model from the viewpoint of its stability.Partially, operational stability deals with the asymptotic nature of the stability notion, whereas biochemists conduct experiments during "finite time." In turn, the stability theory of dynamic systems, so-called Lyapunov stability, offers powerful tools for enhancing biochemical reactions modeling and saving the qualitative behavior of the systems.Moreover, each model is based on a series of assumptions about biochemical interactions, which allows us to check the validity of phenomena that cannot be verified experimentally, for example, mass action law with the distributed or discrete delays, which will be presented later.

Background for lactate measurement
Lactate characteristics: Lactate is an anion of lactic acid and is the final metabolite of the anaerobic breakdown of glucose.It is formed from pyruvate during the processes of glycolysis in the absence of oxygen (Figure 1) (Rabinowitz and Enerbäck, 2020) and is an important substance used in medicine as a marker of hypoxia and a number of other disorders, including diabetes and liver and kidney disorders (Pundir et al., 2016).
Lactate is also actively used in the food industry as an emulsifier, thickener, and acidity regulator.Various salts of lactic acid in the international classification of food additives are numbered E325, E326, E327, E328, and E329.Magnesium lactate is sometimes classified as an antioxidant (Standards, 2021).In addition, lactate can act as an indicator of the activity of bacteria during the fermentation process (Rawoof et al., 2020).It can indicate the freshness and quality of some products-wine (Gamella et al., 2010), juices (Trifirò et al., 1997), etc.
The concentration of lactate in products increases as they spoil due to changes in organoleptic properties.Due to the wide use of lactate as a marker of many processes in medicine and industry, accurate and effective methods of its diagnosis are necessary.

State-of-the-art lactate measurement with the help of biosensors
Biosensors offer effective and sensitive detection methods that can be used in medical institutions to measure the level of lactate in blood, sweat, and other biological fluids.In addition, new methods of lactate detection can be useful for enterprises that must control the processes of manufacturing food products and pharmaceuticals (Rattu et al., 2020).
Classification of lactate biosensors from the viewpoint of a bioselective element: A number of biosensor developments aimed at determining the concentration of lactate are known to be in development.Known biosensors can be divided according to the type of bioselective element: sensors based on lactate oxidase (LOx) and sensors based on lactate dehydrogenase (LHD) (Rassaei et al., 2013).In both cases, the substrate and product of the enzymatic reaction are lactate and pyruvate, respectively.
Basic reactions for lactate biosensors: There is a significant difference in these two reactions.For the LHD reaction, NAD+ is needed as a proton carrier in the reaction of dehydrogenation of lactate to pyruvate (Chatterjee et al., 2023).The lactate oxidase reaction can be much easier because, in it, the role of the proton acceptor is played by oxygen; that is, the use of this enzyme in biosensors does not require additional reagents (Lockridge et al., 1972).
Other variants of biosensors are based on a mixture of LHD and LOx.This configuration of the bioselective element, according to the results obtained by Chaubey et al. (2000), allows measuring lactate at lower concentrations than mono-enzyme biosensors, but this makes the analysis more expensive and is characterized by the complexity of manufacturing.
Classification of lactate biosensors from the viewpoint of the transducer: The main biosensor measurement methods used for lactate analysis are optical-electrochemiluminescence or fluorescence (Pundir and Narwal, 2017) and electrochemical-amperometric, potentiometric, or, less often, conductometric or impedimetric (Rathee et al., 2016).Moreover, various methods of improving the main characteristics of sensors have been applied to known biosensor systems-nanoparticles (Nesakumar et al., 2013;Azzouzi et al., 2015), other nanomaterials (Cui et al., 2007), complex, multi-stage methods of enzyme immobilization (Parra et al., 2006), and multienzyme membranes.
Motivation for lactate biosensor design: Therefore, in connection with the variety of works on the development of biosensor systems for measuring lactate and the prospect of its wide use in various areas of human life, we have concluded that the development of new biosensor methods for determining lactate is an urgent need.
Most existing biosensors are currently not ready for wide implementation and commercialization due to various limitations, such as insufficient sensitivity, selectivity, stability with respect to possible inferents, or too narrow a range of biosensor operation.

Michaelis-Menten model for enzyme kinetics
According to the model, an enzyme E combines with a substrate S to form an enzyme-substrate complex ES, characterized by a rate constant k 1 .The resulting complex can dissociate into E and S (with a rate constant of k −1 ) or transform into a product P with a rate constant of k 2 (Berg et al., 2002).
The speed of the enzyme process is dependent on the ease of formation of the complex of the enzyme with the substrate.For low substrate concentrations, the reaction rate is proportional to the substrate concentration, while at higher concentrations, it tends toward a maximum value and becomes independent of substrate concentration.The general dependence of the rate of an enzyme reaction on substrate concentration is described by an equation called the Michaelis-Menten equation (Eq.1): Here, v is the rate of reaction (mol/s), V max is the maximum reaction rate (mol/s), n S is substrate concentration (mol/dm 3 ), and K M is the Michaelis-Menten constant (mol/dm 3 ).The Michaelis-Menten constant from Eq. 1 is an enzymespecific quantity, dependent on substrate, temperature, and pH and independent of enzyme concentration.This constant is a measure of the affinity of the enzyme for the substrate.The lower the value of the constant, the higher the affinity of the enzyme for the substrate (Berg et al., 2002;Radomska, 2016).

Stability research on enzyme kinetics
Models of enzyme kinetics are based on compartmental systems, which are dynamic systems characterized by a network of interconnected nodes, each representing a reservoir or compartment where resources are stored (Blanchini et al., 2023).The system's behavior is governed by the movement of resources, depicted as flows traveling along the edges connecting these compartments.Compartment-based dynamic systems serve as invaluable models across various disciplines, including physiologically based pharmacokinetics (Martsenyuk et al., 2012;Thompson and Beard, 2012), mathematical epidemiology (Reyné et al., 2022;Martsenyuk et al., 2021a;b), enzyme kinetics (Keener and Sneyd, 2009;Craciun et al., 2020;Martsenyuk et al., 2022), demography (Navarro Valencia et al., 2023), and ecology (Krishna et al., 2024).Stability analysis plays a pivotal role in understanding the behavior of these systems under different conditions (Blanchini et al., 2023).In recent years, researchers have made significant strides in advancing stability research methodologies, particularly in the context of compartmental systems with delay (Martsenyuk et al., 2013;Martsenyuk and Gandzyuk, 2013).
One common approach to stability analysis is linearization, which involves approximating nonlinear systems around equilibrium points.This technique has been extensively utilized in physiologically based pharmacokinetic models to assess the stability of drug distribution processes within the body (Martsenyuk et al., 2012).Lyapunov functions represent another powerful tool for stability analysis, offering a rigorous mathematical framework to prove the stability properties of compartment-based systems.In mathematical epidemiology, Lyapunov functions have been employed to establish global stability of disease-free and endemic equilibria in compartmental models of infectious diseases (Martsenyuk et al., 2021a).

Brief description of the work
Section 2 describes the materials, including the experiment and models used.The methods presented are related to parameter identification and stability research.Section 3 shows the results concerning the parameter identification for models using gammadistributed delay and with two discrete delays.The qualitative analysis includes the existence and positiveness of the solutions, equilibrium states, marginal stability, and numerical research with the help of Poincaré sections.In Section 4, we discuss the results obtained and open problems.
The objective of the work is to offer the flowchart of the lactate biosensor design, including modeling, parameter identification, and stability analysis.

Chemical compounds
The enzyme lactate oxidase obtained from Aerococcus viridans with an activity of 100 units (Sigma, United States) was used to create the biosensor.Bovine serum albumin (fraction V) (BSA) and a 25% aqueous solution of glutaraldehyde (GA), glycerol, and KCl were obtained from Sigma-Aldrich (United States).Stock 500 mM sodium L-lactate solution (Sigma-Aldrich, Switzerland) was used as a substrate.HEPES obtained from Sigma-Aldrich (United States) was used to prepare the buffer solution.Other inorganic compounds used in the work were of domestic production and had a purity level of "h.p." and "p.d.a." Scheme of reactions that underlie the functioning of a biosensor.

Lactate biosensor design
The general scheme of amperometric transducers is shown in Figure 2. Platinum disc electrodes were used as amperometric transducers, manufactured in the laboratory of the Department of Biomolecular Electronics of the Institute of Biomolecular Biology and Geosciences using the following technology: a piece of platinum electrode with a diameter of 0.5 mm and a length of 3 mm was placed in a capillary tube with an outer diameter of 3.5 mm, and then the narrowed end of the capillary tube was sealed in a torch flame.The electrical connection between the platinum and the silver wire conductor was made by low-temperature soldering using Wood's alloy.The open end of the capillary was filled with epoxy resin, with part of the conductor inside the capillary and part outside.A copper contact was soldered to the conductor to connect the transducer to the measuring unit.Amperometric measurements setup is shown in Figure 3.The PalmSens potentiostat (Palm Instruments BV, the Netherlands) was connected to an auxiliary platinum electrode, a silver chloride (Ag/ AgCl) reference electrode, and working electrodes based on platinum disc electrodes.The potentiostat was connected to an 8-channel device (CH-8 multiplexer, Palm Instruments BV, the Netherlands), which allowed it to receive signals simultaneously from several working electrodes or biosensors (up to eight simultaneously).The distance between the auxiliary platinum electrode and all working biosensors during the measurement was the same (approximately 5 mm).
Preparation of bioselective membranes: Bioselective membranes were prepared by immobilizing proteins on the surface of a platinum disc electrode by covalent cross-linking of the enzyme in a bovine serum albumin matrix using glutaraldehyde as a cross-linking agent.The enzyme gel containing 5% lactotoxidase, 5% BSA, and 10%  glycerol was mixed with a 1% glutaraldehyde solution in a 1:1 ratio, after which the mixture was applied to the sensitive area of the platinum disc electrode, and the membrane was air-dried for 20 min.After immobilization, the residual glutaraldehyde and unbound membrane components were washed off the membranes in a buffer solution for 5 min with constant stirring, changing the buffer several times.

Description of amperometric measurements
The measurement was carried out in an open measuring cell with a volume of 2 mL under constant stirring at room temperature.A fixed potential of +0.6 V was applied to the electrodes relative to the chloride silver reference electrode.The working buffer was 25 mM HEPES with pH 7.4.The required concentration of the substrate in the cell was set by adding aliquots of the standard solution, 500 mM sodium L-lactate, to the buffer.
The duration of one response (from the addition of the substrate to the signal output to the baseline) was approximately 4-5 min; between responses, the substrate was washed off the biosensor for 5 min, changing the buffer in the measuring cell several times.All measurements were performed in at least three replicates.

Models used 2.4.1 Model with continuously distributed delays
An application of delayed mass action law to enzyme kinetics was inspired by Brown's model, formulated by Brown (1902), where complex C has a lifetime τ before being decayed.We called the reaction scheme an irreversible one-complex Brown's (IR1CB) mechanism.In Martsenyuk et al. (2022), we offered the following model based on continuously distributed delays: where for confidence level c ∈ (0, 1), we set is the density function of the delay distribution, which was designed in the form of a gamma distribution: where a, m, τ min ≥ 0 are the parameters that determine the corresponding probability density function.Their roles and the ways of estimating were well-studied in Martsenyuk et al. (2022).
The basic idea of Brown's model shown in Model (2) does not include complex C directly but involves the model time τ required for complex forming-destroying.The model parameters were wellstudied by Martsenyuk et al. (2022) and can be used as an initial approximation for the complex-based model in the next subsection.

Models with two discrete delays
The model extends the well-known irreversible one-complex Michaelis-Menten (IR1CMM) mechanism (Keener and Sneyd, 2009) by adding the time durations τ 1 and τ 2 required for the entire fulfillment of the forward reactions.Mathematically, it corresponds to the time delays within dynamic systems.Hence, we consider the following set of elementary reactions: E + P, which we call the irreversible one-complex with two delays Michaelis-Menthen (IR1C2DMM) mechanism.Based on the general approach described by Craciun et al. ( 2020), it yields the delayed model For the solutions of Eq. 4, elements of which are the vector functions n S , n E , n C , n P ∈ C 1 ([−τ max ], 0], R 4 ), we consider the following initial conditions: (5) The value of n P (t) can be found by direct integration, namely, 2.5 Methods

Parameter identification
As a result of the amperometric measurements described in Section 2, we obtain responses in the form of current I(t).We propose to use the relation of the current I(t) with n P (t) as In Martsenyuk et al. (2022), this relationship was evidenced for the specific conductance κ(t) of conductometric biosensors.Mathematical modeling of conductometric biosensors in terms of conductivity is presented in detail by Zouaoui et al. (2022).Provided fixed potential, we assume the linear dependence of the specific conductance and the current.So, we follow Eq. 7. Numerical modeling regarding amperometric biosensors is displayed by Simelevicius et al. (2012) andHashem Zadeh et al. (2020).
In the following, we will denote product concentrations obtained as the solutions of the models as n P,pred (t).In turn, the corresponding values of the current due to Eq. 7 will be I pred (t).On the other hand, let I exp (t) be the values of responses received as a result of the experiments.
The proposed parameter identification uses the currents I exp,j (t i ) and I pred,j (t i ), j 1, m, experimentally and numerically obtained at the time instances t i , i 1, N1 for given initial substrate concentrations n S (0) = n S,j .
Our goal is to estimate the parameters when applying Model (2) or (4), respectively.We will estimate the parameters Π of the models solving the problems of optimization: Here, is the target function, and are inequality constraints, where Θ is a null vector, and Π lower and Π upper are lower and upper bounds for the parameter values of corresponding dimensions.
An algorithm for the estimation of Π IR1CB was described in detail by Martsenyuk et al. (2022).Here, we also apply it for the estimation of Π IR1C2DMM .Moreover, the values of a, m, and τ min, obtained as a result of the parameter identification in Eq, 2, will be used in order to set the initial estimate for τ 1 and τ 2 when estimating Eq. 4. Note that the mean value of time delay for Model (2) can be calculated as So, we set the initial values of delays from Eq. 4 such that τ 1 + τ 2 ≈ E(τ).

Stability research using linearization technique
We will conduct research on local stability using the linearization technique.In this case, stability conditions are constructed based on a characteristic quasi-polynomial.The signs of the real parts of its roots are decisive for making conclusions about stability.
Namely, if all roots lie on the open left-half plane, then the equilibrium is locally asymptotically stable.If some of the roots have positive real parts, then the equilibrium is unstable.We focus our attention on the case when the roots lie on a closed left-half plane; that is, we have some simple roots lying on an imaginary axis.Such a situation is known as marginal stability.We distinguish this situation from the case of the pair of purely imaginary roots corresponding to periodic solutions.
When investigating the characteristic quasi-polynomial, a Padé approximant will be used in the form allowing us to approximate the characteristic quasi-polynomial with the help of rational functions (Baker and Graves-Morris, 1996).

Bifurcation plots using Poincaré section
We use the Poincaré section technique to study the qualitative behavior of the models developed, which was primarily applied to the compartmental model by Martsenyuk et al. (2021a).
To begin, we gain a thorough understanding of Model (4), including its equations and parameters.We select parameters τ 1 and τ 2 to vary and the variables n S , n E , and n C to focus on.Then, we simulate the system across different parameter values, generating time series data for the chosen variables.
Poincaré sections are constructed by intersecting the trajectory of the system with a defined plane n E = d in the phase space, where d = (min t n E (t) + max t n E (t))/2.This section is determined by specific criteria; namely, we choose such points (n + S , d, n + C ) such that crossing a plane n E = d will happen at sequential time instances t + + T, t + + 2T, t + + 3T, . . ., where T is a period value; that is, n + S n S (t + ) n S (t + + T) n S (t + + 2T) . . . .We plot the sampled points in an (n S , n E )-space to visualize the Poincaré sections.
Repeating this process for various parameter values (τ 1 , τ 2 ), we observe how the Poincaré sections change.Analyzing the patterns in the sections, we understand the system's behavior, including the emergence of periodic orbits, chaotic dynamics, or transitions between different states.
Finally, we summarize the results by constructing a bifurcation plot combined with the corresponding n S , n E and (n S , n P ) phase plots, showing how features of the Poincaré sections vary with parameter values.This comprehensive approach allows us to systematically explore the system's behavior and gain insights into its dynamics.

Parameter identification for a model using a gamma-distributed delay
The parameter identification technique mentioned in Section 2.5.1 was used for Model (2).The training data correspond to a set of time series of currents corresponding to given initial substrate concentrations n S (0) equal to 0.1 mM, 0.5 mM, 1.0 mM, and 2.5 mM sequentially.
The goal is to estimate the parameters The initial values for the estimation were chosen as Π init (see Table 1).The lower and upper bounds for Π IR1CB are shown in Table 1 as Π lower and Π upper , respectively.
The most valuable information obtained from Model (2) is the density of the distributed delay distribution (see Figure 4), which will be used in further estimations.
The comparison of predicted and expected values of the currents for the optimal set of the parameters Π opt is shown in Figure 5.We see that the parameter values, being optimal concerning the cost criterion (8), enable us to find the solution of Model (2) closest to the expected currents for the initial substrate concentration of 1.0 mM.For the smaller initial values of the substrate, we have predicted values smaller than the expected ones, whereas for those bigger than 1.0 mM, the predicted values are larger than the experimental ones.The explanation of such an effect lies in the special kind of stability of the model known as marginal stability, which will be evidenced further.

Parameter identification for a model with two discrete delays
Model (4) requires the estimating the parameters They were obtained as a result of the solution of the optimization problem shown in Eqs 8-10.The training data described in Section 3.1 were used.
The initial values for the estimation were chosen as Π init (see Table 2).The lower and upper bounds for Π IR1C2DMM are shown in Table 2 as Π lower and Π upper , respectively.
The comparison of predicted and expected values of the currents for the optimal set of the parameters Π opt for Model (4) is shown in Figure 6.We see similar tendencies of expected and predicted values as we did for Model (2), with gammadistributed delays.
Henceforth, we will focus our attention on the positiveness of the solution of the system shown in Eq. 4.
The positiveness of n P (t) follows directly from the positiveness of n C (t).So, we prove the positiveness of n S (t), n E (t), and n C (t) by contradiction.

Case without delays
First, we demonstrate the positiveness for Model (4) without delays, that is, if Let us assume, for the sake of contradiction, that there is the smallest value among t c , t e , and t s delivering non-positive solutions.Consider them sequentially.
Let t c > 0 be the smallest instance of time such that n C (t c ) = 0. From the first and second lines of Eq. 12, we get It implies that Density function given in 3 for the distributed delay τ value from Model Eq. ( 2).
From the third line of Eq. 12, we have Hence, By continuing Eq. 13, we have n C (t c ) > 0, which contradicts the initial assumption.
Let t e be the smallest instant that n E (t e ) = 0. From the third part of Eq. 12, it follows that and n C (t) > 0 for t ∈ [0, t e ).
In turn, from the second part of Eq. 12, we have and Let t s be the smallest instant that n S (t s ) = 0. From the second line of Eq. 12, we get that n C (t) > 0, t ∈ [0, t s ).Furthermore, from the first equation, we have and Comparison of the predicted and expected currents for Model (2) at the optimal values of the parameters Π opt corresponding to given initial substrate concentrations n S (0): (A) 0.1 mM, (B) 0.5 mM, (C) 1.0 mM, and (D) 2.5 mM.
For the reasons given, we see that the solution of Eq. 12 exists and is positive for any positive initial values (n S (0), n E (0), n C (0), n P (0)) > 0.

Two discrete delays
It is natural to assume that the solutions of Eq. 4 are still positive for some sufficiently small τ 1 and τ 2 .We aim to offer the conditions of positiveness.
The conditions will be based on the notion of the delayed exponential function.Given the value x ∈ R and delay τ > 0, the delayed exponential function is called (Angstmann et al., 2023) where Θ(•) is the Heaviside function.This function has the property that d dx e λx τ λe λ(x−τ) τ .Note that contrary to the "undelayed" exponential function, e x τ can also accept negative values.To apply the approach of contradictions shown above, let t c , t e , and t s be instances delivering non-positive solutions.
If t c > 0 be the smallest instance that n C (t c ) = 0, then from the first and second lines of Eq. 4 Hence, , where n max S ≔ max u n S (u).From the third line of Eq. 4, we have Consider t = t c .Hence, and we see that Following the proof in the previous case, we conclude that the solutions of Eq. 4 are positive if rate parameters, initial conditions, and delays τ 1 , τ 2 are such that for any t > 0 e −k1nE 0 t τ1 > 0, e −k1n max S t τ1 > 0, e −k2nE 0 t τ2 > 0.

Equilibrium of the system
Let ( n S , n E , n C , n P ) be the equilibrium of Model (4).It should satisfy When analyzing Eq. 4, we see that there is a conserved quantity of enzyme because d(nE+nC) dt ≡ 0. Hence, n E (t) + n C (t) ≡ n E0 , where n E0 is the total amount of available enzyme.Thus, from the last equality of Section 3.3.2.1, we conclude that n E n E0 and n S 0. So, we have the unique equilibrium of Model (4), which is the substrate and complex free equilibrium (SCFE).
n S 0, n E n E0 , n C 0, n P const.
The value of n P is bounded and can be determined by Eq. 6.Note that it is related to the initial values of both n S and n E .n P is not included in the right-hand sides of Eq. 4.

Marginal stability
The Jacobian matrix at SCFE for Model (4) is given by In Figure 9, we depict the bifurcation diagram for the variable n S resulting from the parameter variation of τ 1 .It reveals that the equilibrium state behaves as a stable node within the range of τ 1 values from 0 to 25. Subsequently, a bifurcation occurs, leading to a transition to a limit cycle.This period-doubling phenomenon is further supported by Figure 10, which displays the (n S , n E )-phase plane for different initial states.
For τ 1 values exceeding 32, the system's behavior appears likely to become chaotic.This transition is occurring through period doubling, although further investigation is warranted.It is noteworthy that the system exhibits marginal stability, as evident from Figure 11.This figure illustrates how different initial values of n S (0) lead the solution n P to converge toward either a stable node or a limited cycle.
Special attention should be paid to the positivity of the solutions, as we observed from the numerical modeling.

Discussion
The use of gamma-distributed delays and discretely distributed delay models is simultaneously important because the parameter values from the IR1CB model can be used for the IR1C2DMM model.The gamma-distributed delays model allows us to estimate the mean value of the delay easily.Moreover, based on the results of numerical experiments, we can conclude that for the optimal values of the parameters Π opt (Tables 1, 2), we have that τ 1 + τ 2 ≈ E(τ), where τ 1 , τ 2 are the discrete delays from Eq. 4, and τ is the distributed delay from Eq. 2 (see also Figure 4 for delay density distribution).
The use of the Levenberg-Marquardt algorithm is essentially determined by choosing the initial approximation of the parameter values.It can be improved by applying AI techniques and constructing and tuning neural networks of the appropriate architecture for future consideration.
Artificial intelligence (AI) models like recurrent NNs can also be used to model enzyme-substrate interaction (we will try this in the future).The advantage of an AI model is that it makes more accurate predictions.The disadvantages of using AI are (1) the black box problem, meaning that the models are not based on any biochemical assumptions but only neural network architecture and (2) the overfitting problem, meaning that when we try to fit the outputs as accurately as we can, we may not see the general tendencies in reactions.
The system is characterized by marginal stability, which is between Lyapunov stability and instability.It corresponds to the objective of an electrochemical biosensor as a measuring device.Each initial state, including substrate concentration changes, has its "own" concentration of the product at an infinite time.We can reformulate it as the definition of operating stability.The question arises of how to apply and interpret the marginal stability condition.As was shown, the marginal stability condition can be reduced by checking the real parts of the polynomial roots.We see that the roots depend on the reaction rate parameters (k 1 , k 2 , k −1 ), time delays τ 1 , τ 2 , and the initial concentration of n E0 .

Conclusion
The equilibrium of the system was demonstrated in this article.From the viewpoint of enzymatic reactions and based on the proposed model, we saw that the system would be in equilibrium in the case of the absence of substrate and complex, that is, n S ≡ 0, n C ≡ 0, as well as the concentration of the enzyme equal to the total amount of available enzyme.The concentration of the product should be a constant value calculated in accordance with Eq. 6 as the area under curve n C (t), that is, n P ≡ const.In addition, n P is determined by the initial value n S (0).
In the work, we also considered a model with two discrete delays where τ 1 and τ 2 are time durations required for the entire fulfillment of forwarded reactions.To be precise, τ 1 is the time needed for the substrate to bind with the enzyme (time needed for complex formation), and τ 2 is the time needed for the complex to break down into enzyme and product.It was numerically shown that increasing those two parameters could result in a loss of stability.Suppose we have optimal values τ 0 1 and τ 0 2 and they correspond to a stable solution.If τ 1 > τ 0 1 , more time is required for forming the complex at S + E → k1,τ1 C. In other words, less complex will be produced compared with the optimal value.In turn, if τ 2 > τ 0 2 , this means that more time will be needed for the C → k1,τ2 E + P reaction, leading to less enzyme and product than in the steady state.We remember that in the steady state, it has been proven that we have no complex, and the enzyme content is maximal.On the other hand, an increase in τ 2 results in having more complex but less product.This will affect the final content of the product, that is, the biosensor index.
It can also be noted that no unique parameters are set for fitting the model to the expected data at different values (Figures 5, 6), as the rates are not constant but can be functions of the substance concentrations.In particular, the times required for the reactions S + E → k1,τ1 C and C → k1,τ2 E + P are not constant.They likely depend on the concentrations of S, E (for τ 1 ), and C (for τ 2 ).These data should be determined experimentally.Moreover, the models can be extended by involving other complexes created by reactions with other biosensor auxiliary components.For example, BSA is known to create a complex with an enzyme known as CLEA (Shah et al., 2006).In other words, incorporating additional variables and increasing the system dimension can improve the model fitting.
These findings provide valuable insights into the system's dynamics and highlight the importance of parameter exploration in understanding the qualitative behavior of enzyme kinetics models.

FIGURE 2
FIGURE 2 General scheme of amperometric transducers.(A) Schematic view of the amperometric transducer.(B) Photograph of an amperometric transducer.(C) Scheme of the sensitive area of the transducer.(D) Photograph of the sensitive area of a transducer.1-sensitive area, 2 -enzyme membrane, 3-platinum wire, 4-inner conductor (silver wire), electrical connection using 5-low-melting Wood's alloy, 6-epoxy resin, and 7-contact area.

FIGURE 6
FIGURE 6Comparison of the predicted and expected currents for Model (4) at the optimal values of the parameters Π opt corresponding to the given initial substrate concentrations n S (0): (A) 0.1 mM, (B) 0.5 mM, (C) 1.0 mM, and (D) 2.5 mM.